%% Stock Return Extrapolation, Option Prices, and Variance Risk Premium
% By Adem Atmaz
% This file generates Table 8 Panel B: Leverage effect and the variance risk premium
% It also requires 1 function file: olshac.m
% Typically it takes less than 2 minutes to run

clear all;  clc;
format compact
format

% Parameter Values follow from  "Atmaz_Table_4_ParameterValues.m"
alpha=0.50;
theta=0.25;
mu=0.0212;
kappa=0.3567;
zeta=0.0044;
sigma=0.1625;
rho=0;   % vary correlation to rho=0.415
r=0.0056;

D0=100;
a=2.75;
vec_theta=[0, 0.25, 0.50, 0.75];
N_th=length(vec_theta);



% Simulations and Regressions
tau=1/12;  % variance horizon in years
h=3;       % predictive regression holding return horizon in months
N_sim=1000;                       % number of simulations We use 1000
S_end_year=27;                    % Match sample size
int = 1/252;                      % time interval. Use 1/(k*252) for k points in each day in  simulation
time_grid = 0:int:S_end_year;     % Time vector from 0 to end year
N_time=length(time_grid);
N_rets=N_time-1;                % Number of returns in each path
H=1/12;     % Data frequency in the regression e.g. 1/12 for monthly
N_pp=H/int; % number of data points for each data point in regression: 21 for month

% omitting the first burn out period
keep_final_years=S_end_year; % keep all of it
Obs=keep_final_years*12;

% 1 path vectors
vec_D = zeros(1,N_time); % Allocate output vector
vec_V = zeros(1,N_time); % Allocate output vector
vec_X = zeros(1,N_time); % Allocate output vector

% All paths monthly matrices
N_month=(12*S_end_year)+1;

SM_V = zeros(1,N_month); % Allocate output vector
SM_X = zeros(1,N_month); % Allocate output vector
SM_D = zeros(1,N_month); % Allocate output vector
SM_S = zeros(1,N_month); % Allocate output vector
SM_VSR = zeros(1,N_month); % Allocate output vector
SM_VRP = zeros(1,N_month); % Allocate output vector
SM_RV  = zeros(1,N_month); % Allocate output vector

SM_Cov= zeros(1,N_month-1); % Allocate output vector
SM_Corr= zeros(1,N_month-1); % Allocate output vector

% HIGH AND LOW Vectors
vec_beta_High=      zeros(N_sim,N_th);
vec_tstat_hac_High= zeros(N_sim,N_th);
vec_R2_High=    zeros(N_sim,N_th);
vec_beta_Low=       zeros(N_sim,N_th);
vec_tstat_hac_Low=  zeros(N_sim,N_th);
vec_R2_Low=     zeros(N_sim,N_th);
vec_High_corr=      zeros(N_sim,N_th);
vec_Low_corr=       zeros(N_sim,N_th);

% Begin simulation and regression
rng('default'); % fix the random number generator  seed  for replication
vec_Z1=randn(N_sim,N_time);
vec_Z2=randn(N_sim,N_time);
tic
for s = 1:N_sim
    
    for th=1:N_th
        theta=vec_theta(th);
        % Stock price quantities
        M_DELTA=((kappa+(rho*sigma)*(1+theta))^2)-((2+theta)*(1+theta)*sigma^2);
        M_c=theta/(alpha*(1+theta));
        M_b=-(2+theta)/(kappa+((rho*sigma)*(1+theta))+sqrt(M_DELTA));
        M_a=a;
        M_sigma_1=(1+(rho*sigma*M_b))/(1-alpha*M_c);
        M_sigma_2=(sqrt(1-rho^2)*sigma*M_b)/(1-alpha*M_c);
        M_Var_con=(M_sigma_1^2)+(M_sigma_2^2);
        M_Cov_con=(rho+(sigma*M_b))/(1-alpha*M_c);
        
        % VSR and VRP quantities
        M_zeta_v=zeta*M_Var_con;
        M_kappa_v=kappa+(sigma*M_Cov_con);
        M_eta_v=(theta/M_b)*(1-(alpha*M_c))*M_Var_con;
        
        v0=M_zeta_v/kappa;
        M_RN_LRM_vt=(M_zeta_v+(r*M_eta_v))/(M_kappa_v+(0.5*M_eta_v));
        
        tt1=M_kappa_v+alpha;
        kappa1=(0.5*tt1)-(0.5*sqrt((tt1^2)-(4*alpha*(M_kappa_v+(0.5*M_eta_v)))));
        kappa2=(0.5*tt1)+(0.5*sqrt((tt1^2)-(4*alpha*(M_kappa_v+(0.5*M_eta_v)))));
        EE0=(1-exp(-kappa*tau))/(kappa*tau);
        EE1=(1-exp(-kappa1*tau))/(kappa1*tau);
        EE2=(1-exp(-kappa2*tau))/(kappa2*tau);
        
        % RV A and B
        RV_A=v0*(1-EE0);
        RV_B=EE0;
        
        % VSR A and B and C
        VSR_A=(M_RN_LRM_vt*(1-EE1-((kappa1/(kappa2-kappa1))*(EE1-EE2))))...
            +((M_zeta_v/(kappa2-kappa1))*(EE1-EE2));
        VSR_B=EE1-(((M_kappa_v-kappa1)/(kappa2-kappa1))*(EE1-EE2));
        VSR_C=(M_eta_v/(kappa2-kappa1))*(EE1-EE2);
        
        % VRP A and B and C
        VRP_A=(v0*(1-EE0))-VSR_A;
        VRP_B=EE0-VSR_B;
        VRP_C=-VSR_C;
        
        % Initalization of each path
        D0=100;
        V0=zeta/kappa; % long run mean of  fundamental  variance
        X0=(r+(0.5*M_Var_con*V0))/(1+theta); %  inital  X is equal to its long run mean
        
        vec_D(1)=D0;
        vec_V(1)=V0;
        vec_X(1)=X0;
        
        for j = 1:N_rets
            Z1=vec_Z1(s,j);
            Z2=vec_Z2(s,j);
            vec_V(j+1) = max(0, vec_V(j))+((zeta-(kappa*max(0, vec_V(j))))*int)...
                +(sigma*rho*sqrt(max(0, vec_V(j)))*sqrt(int)*Z1)...
                +(sigma*sqrt(1-rho^2)*sqrt(max(0, vec_V(j)))*sqrt(int)*Z2);
            vec_X(j+1) = vec_X(j)+((r+(0.5*M_Var_con*vec_V(j))-(1+theta)*vec_X(j))*alpha*int)...
                +(alpha*M_sigma_1*sqrt(max(0, vec_V(j)))*sqrt(int)*Z1)...
                +(alpha*M_sigma_2*sqrt(max(0, vec_V(j)))*sqrt(int)*Z2);
            vec_D(j+1) = vec_D(j)+((mu*vec_D(j))*int)+(vec_D(j)*sqrt(max(0, vec_V(j)))*sqrt(int)*Z1);
            
        end
        
        %  get daily returns and variance changes
        SM_daily_S=vec_D.*exp(M_a+(M_b*vec_V)+(M_c*vec_X));
        SM_daily_Rets=(SM_daily_S(2:end)./SM_daily_S(1:end-1))-1;
        SM_daily_dvar=M_Var_con*(vec_V(2:end)-vec_V(1:end-1));
        
        % Get monthly conditional correlations from daily values
        for j = 1:N_month-1
            j_str=(N_pp*(j-1))+1;
            j_end=j_str+N_pp-1;
            aa=cov(SM_daily_dvar(j_str:j_end),SM_daily_Rets(j_str:j_end));
            SM_Cov(j)=aa(1,2)*10000; % interms of percentage squared
            SM_Corr(j)=corr(SM_daily_dvar(j_str:j_end)',SM_daily_Rets(j_str:j_end)');
            
        end
        LT=length(SM_Corr);
        
        % Get end of month values for each process from daily values
        SM_V= vec_V(1:N_pp:end);
        SM_X= vec_X(1:N_pp:end);
        SM_D= vec_D(1:N_pp:end);
        SM_S= SM_D.*exp(M_a+(M_b*SM_V)+(M_c*SM_X)); %stock price
        
        
        % +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        %  #######################      instantenous  RV       ################
        %  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        RV_A= M_zeta_v;
        RV_B=1-kappa;
        
        % +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        %  #######################   instantenous     VSR       ###############
        %  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        VSR_A=M_zeta_v;
        VSR_B=1-M_kappa_v;
        VSR_C=M_eta_v;
        
        % simulated VSR VRP and RV
        SM_RV=(RV_A+(RV_B*M_Var_con.*SM_V))*(10000/12); %normalized to get monthly values in percentage squared
        SM_VSR=(VSR_A+(VSR_B*M_Var_con.*SM_V)+(VSR_C.*SM_X))*(10000/12);  %normalized to get monthly values in percentage squared
        SM_VRP=SM_RV-SM_VSR;        
        
        % Keep only the last observations
        SM_RV=SM_RV(end-Obs:end);
        SM_VRP=SM_VRP(end-Obs:end);
        SM_S=SM_S(end-Obs:end);
        SM_D=SM_D(end-Obs:end);        
      
        % get h period ahead log returns
        SM_Ret_h=log(SM_S(1+h:end))-log(SM_S(1:end-h)); % log ret
        %SM_Ret_h=(SM_S(1+h:end)./SM_S(1:end-h))-1;      % returns
        SM_Ret_h_ann=(12/h)*SM_Ret_h*100;  %annualized in percentage terms
        
        High_val_corr=quantile(SM_Corr,0.33);
        Low_val_corr=quantile(SM_Corr,0.67);       
        
        High_corr_=zeros(1,LT);
        Low_corr_=zeros(1,LT);
        High_VRP_=zeros(1,LT);
        High_Ret_=zeros(1,LT);
        Low_VRP_=zeros(1,LT);
        Low_Ret_=zeros(1,LT);
        
        for j = 1:LT-h
            corr_j=SM_Corr(j);
            if  corr_j<High_val_corr   % get high correlation vectors
                High_corr_(j)=corr_j;
                High_VRP_(j)=SM_VRP(j+1);
                High_Ret_(j)=SM_Ret_h_ann(j+1);
            end
            
            if  corr_j>Low_val_corr   % get low correlation vectors
                Low_corr_(j)=corr_j;
                Low_VRP_(j)=SM_VRP(j+1);
                Low_Ret_(j)=SM_Ret_h_ann(j+1);
            end
            
        end % end j for each month
        
        % get nonzero values
        vec_High_corr(s,th)=mean(High_corr_(High_corr_~=0));
        vec_Low_corr(s,th) =mean(Low_corr_(Low_corr_~=0));
        High_VRP = High_VRP_(High_VRP_~=0);
        High_Ret = High_Ret_(High_Ret_~=0);
        Low_VRP  = Low_VRP_(Low_VRP_~=0);
        Low_Ret  = Low_Ret_(Low_Ret_~=0);        
        LT_High=length(High_VRP);
        LT_Low=length(Low_VRP);        
        
        % Run future returns on VRP for high correlation months
        reg_Y=High_Ret';
        reg_X=[ones(LT_High,1) High_VRP'];
        lag=round(0.75*(LT_High^(1/3)));
        
        [beta, R2, R2adj, X2_NW, X2_HH, X2_R, std_NW, std_HH, std_R, t_NW, t_HH, t_R]=olshac(reg_Y,reg_X,lag,lag);
        
        vec_beta_High(s,th)=beta(2);
        vec_tstat_hac_High(s,th)=t_NW(2);
        vec_R2_High(s,th) = R2;        
        
        % Run future returns on VRP for Low correlation months
        reg_Y=Low_Ret';
        reg_X=[ones(LT_Low,1) Low_VRP'];
        lag=round(0.75*(LT_Low^(1/3)));
        
        [beta, R2, R2adj, X2_NW, X2_HH, X2_R, std_NW, std_HH, std_R, t_NW, t_HH, t_R]=olshac(reg_Y,reg_X,lag,lag);
        
        vec_beta_Low(s,th)=beta(2);
        vec_tstat_hac_Low(s,th)=t_NW(2);
        vec_R2_Low(s,th) = R2;
        
        
    end   % end theta
    
end    % end sim
clc;
toc

High_R2=round(mean(vec_R2_High)*100,2)
Low_R2=round(mean(vec_R2_Low)*100,2)

clear alpha D0 EE0 EE1 EE2 EEv j jj kappa kappa1 kappa2 m M_a M_b M_c;
clear M_corr_sv M_Cov_con M_DELTA M_eta_v M_kappa_v M_Ob_LRM_vt M_r M_RN_LRM_vt;
clear M_sigma_1 M_sigma_2 M_sigma_v  M_Var_con M_zeta_v mu n  sigma;
clear theta tt1 v V0 Z1 Z2 zeta;
clear vec_Z1 vec_Z2;

%% THE END










